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IV 


ABSTRACT 


This  paper  studies  alternate  ways  to  model  density  diseontinuity  in  split-step  Fourier 
parabolie  equation  models.  The  Monterey-Miami  Parabolic  Equation  model  is  used  to 
implement  an  alternative  to  the  effective  index  term  in  the  smoothing  function  and  a  split- 
step  Fourier/Finite  Difference  hybrid  technique.  The  model  is  shown  to  converge  to  a 
stable  solution  that  is  slightly  lower  than  the  benchmark  solution.  A  range  step  size  of 
approximately  one  wavelength  is  shown  to  provide  the  closest  approximation  to  the 
benchmark  solution.  Acceptable  solutions  are  obtained  with  large  depth  grid  sizes  for  the 
alternate  smoothing  function.  Smaller  depth  grid  sizes  are  necessary  for  accurate 
solutions  when  using  the  hybrid  implementation  technique.  The  effect  of  reference  sound 
speed  is  shown  to  minimize  the  phase  error  present  when  the  models  are  used  in  the 
presence  of  a  strong  density  discontinuity. 
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I.  INTRODUCTION 


A,  BACKGROUND 

Parabolic  equation  (PE)  methods  have  been  used  in  the  field  of  physios  for  many 
years  to  reduoe  the  two-way,  elliptic  wave  equation  to  a  simpler,  one-way  equation. 
Although  originally  utilized  in  the  field  of  plasma  physios,  Hardin  and  Tapper!  [1]  were 
the  first  to  introduce  the  PE  method  for  solving  underwater  acoustic  propagation  in  the 
early  1970s.  The  PE  approximation  to  the  wave  equation  is  first  order  differential  in 
range  (r),  and  therefore  was  able  to  be  numerically  solved  by  using  a  one-way  marohing 
algorithm,  such  as  the  split-step  Eourier  (SSE)  algorithm.  The  SSE  teohnique  works  by 
first  sampling  the  acoustic  field  on  a  depth  grid,  i.e.,  with  uniform  sampling  Az,  at  a 
single  range  r .  The  depth  grid  then  undergoes  two  fast  Eourier  transforms  (EETs)  and 
two  veetor  multiplieations  in  order  to  obtain  the  solution  at  the  next  range  step  r  -I-  Ar . 
The  process  is  repeated  at  eaeh  range  step  until  the  entire  acoustie  field  is  modeled  [2]. 
The  major  advantages  of  the  SSE  technique  are  its  computational  efficiency,  stability,  and 
simplicity. 

Alternately,  finite-differenee  (ED)  techniques  may  be  used  to  solve  PEs. 
Numerieal  implementation  of  the  ED  algorithm  is  similar  to  SSE.  The  aeoustic  field  is 
first  specified  at  all  depths  at  an  initial  range.  A  finite-differenee  approximation  based  on 
a  Taylor  series  expansion  of  the  derivative  terms  is  then  directly  applied  to  the  density- 
dependent,  transverse  differential  operator.  The  solution  is  obtained  by  solving  a 
tridiagonal  system  of  equations  with  dimensions  equal  to  the  number  of  depth  samples  at 
eaeh  range  step  [3].  By  directly  applying  the  finite-differenee  approximation  to  eaeh 
element  in  the  depth  mesh,  the  ED  teehnique  avoids  having  to  use  a  smoothing  function 
and  can  be  used  in  a  variety  of  ocean  environments  but  at  the  eost  of  inereased 
computational  complexity. 

A  third  method  of  solving  PEs  is  the  finite-element  (EE)  method.  This  method 
discretizes  the  oeean  environment  into  multiple  nodes  using  Galerkin’s  method  [4].  Then 
by  applying  continuity  of  pressure  and  the  normal  eomponent  of  particle  veloeity  to  the 
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boundary  conditions  of  each  node,  a  tridiagonal  system  of  equations  is  established  that 
ean  be  solved  through  well-defined  methods.  By  discretizing  the  environment,  the  FE 
method  has  the  ability  to  increase  or  deerease  the  number  of  nodes  within  a  loeal  depth 
range,  which  leads  to  more  detailed  or  faster  ealeulations  than  possible  by  using  FD 
teehniques  [5].  However,  due  to  the  eomplexity  of  FE  sehemes  and  their  assoeiated 
stability  issues,  FD  and  SSF  teehniques  are  the  preferred  method  in  the  modeling 
eommunity  [6]. 

The  1981  PE  Workshop  [7]  proved  Harding  and  Tappert’s  SSF  technique  was 
able  to  efficiently  and  aeeurately  model  the  aeoustie  field  in  long-range,  narrow-angle 
problems  with  negligible  bottom  interaetions.  More  advaneed  SSF  methods  based  on 
higher-order  operator  splittings  of  the  so-ealled  “square  root  operator,”  sueh  as 
introdueed  by  Thomson  and  Chapman  [8],  also  greatly  improved  the  accuracy  of  the 
approaeh  for  wide-angle  solutions. 

In  areas  with  strong  bottom  interaetions,  sueh  as  in  shallow  water,  the  SSF 
teehnique  has  been  shown  to  be  less  aeeurate.  This  is  due  to  the  fact  that  the  propagator 
term  in  the  SSF  algorithm  eontains  density  derivatives  in  the  effeetive  index  of  refraetion 
term,  which  require  the  use  of  a  smoothing  function  to  approximate  the  density  profile  in 
the  vieinity  of  the  oeean  bottom  [9].  This  aecuraey  of  the  smoothing  funetion’s 
approximation  is  eonsidered  to  be  the  major  weakness  to  SSF  based  models. 

In  1997,  Yeviek  and  Thomson  [10]  proposed  a  hybrid  SSF/FD  method.  Their 
initial  motivation  was  to  introduee  a  small  set  of  equations  (generally  five)  that  treat  the 
density  discontinuity  directly  through  a  finite  differenee  approximation  rather  than 
relying  on  the  mixing  funetion  previously  employed,  thereby  improving  the  aecuraey  of 
the  solution.  The  remainder  of  the  ealeulations  and  marching  algorithm  are  performed  via 
SSF  teehniques.  This  method  eapitalizes  on  the  aeeuraey  of  the  FD  teehnique  while 
utilizing  the  effieieney  of  a  SSF  marehing  algorithm.  This  method  has  shown  aeeurate 
results  and  is  effieient  enough  for  use  in  two  and  three  dimensions. 

In  addition,  they  were  able  to  show  some  improved  phase  aeeuraey  (through  less 
sensitivity  to  the  referenee  sound  speed)  by  incorporating  additional  higher  order  terms  in 
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the  operator  splitting  of  the  square  root  operator.  These  additional  terms  were  also  shown 
to  lend  themselves  to  finite  differenee  solutions.  In  this  ease,  however,  the  applieation  of 
the  set  of  finite  differenee  equations  extended  over  the  entire  depth  eomputational 
domain. 

B,  PROBLEM  STATEMENT 

The  Monterey-Miami  Parabolie  Equation  (MMPE)  [11]  method  utilizes  a  SSE 
teehnique  to  model  the  propagation  of  the  aeoustic  field  underwater.  Versions  of  MMPE 
have  been  developed  for  both  2-D  and  3-D,  single  frequeney  and  broadband  ealculations, 
and  ineluding  various  types  of  environmental  fluctuations,  such  as  internal  waves,  surface 
roughness,  3-D  bathymetric  variations,  turbulence,  and  others.  The  MMPE  model  utilizes 
the  wide-angle  PE  (WAPE)  approximation  of  Thomson  and  Chapman  [8],  and  deals  with 
the  density  discontinuity  by  using  a  hyperbolic  tangent  mixing  function  approximation  to 
the  Heaviside  step  function  at  the  bottom  interface. 

This  model  has  been  shown  to  be  accurate  in  deep-ocean  problems  [12]. 
However,  MMPE’s  treatment  of  the  bottom  condition  can  become  problematic  in  shallow 
water  environments  where  bottom  interactions  play  a  significant  role.  In  this  thesis,  we 
will  consider  solutions  to  the  acoustic  propagation  in  a  simple,  Pekeris  waveguide,  where 
the  environment  is  defined  as  having  a  lossless  water  column  with  constant  sound  speed 
of  1500  m/s  and  a  bottom  depth  of  300  m,  below  which  is  a  homogeneous  half-space  of 
sound  speed  1700  m/s,  density  1.5  g/cm3,  and  attenuation  0.5  dB/A.  The  source  is 
modeled  as  an  omnidirectional  point  source  at  depth  180  m  transmitting  a  single¬ 
frequency,  continuous  wave  (CW)  signal  at  100  Hz.  Eor  comparison  purposes,  a 
transmission  loss  trace  will  be  extracted  at  100  m  depth  out  to  a  maximum  range  of 
20  km. 

When  compared  to  a  coupled-mode  model,  known  as  Couple  (and  considered 
capable  of  generating  benchmark  quality  solutions)  and  a  finite-element  parabolic 
equation  (EEPE)  model,  the  MMPE  model  is  quite  accurate  at  shorter  ranges,  as  can  be 
seen  in  Eigure  1.  However,  at  about  5  km  a  phase  error  begins  to  appear  and  accumulate 
as  the  range  is  increased.  The  objective  of  this  work  is  to  study  multiple  implementation 
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options  for  modeling  the  density  diseontinuity  in  MMPE  in  order  to  improve  the  model’s 
long-range  aceuraey.  Insights  into  the  eause  of  this  phase  error  aeeumulation  with  range 
will  also  be  explored. 


Figure  1.  Plot  of  TL  (dB)  vs.  Range  for  Couple,  FEPE,  and  MMPE. 
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II.  THEORY 


A,  PARABOLIC  EQUATION  DERIVATION 

We  begin  with  the  Helmholtz  equation  for  aeoustic  pressure,  in 

eylindrieal  coordinates  for  a  time-dependent  harmonic  field, 


lA 

r  dr 


dp 

dr 


+  p(z)- 
^  ’  dz 


1  dp 
p{z)  dz 


+  kln^{r,z)p  =  0 


(1) 


where  Atq  =  —  is  the  reference  wavenumber,  n{r,  z)  =  — - —  is  the  acoustic  index  of 
Co  c(r,z) 

refraction,  Cq  is  the  reference  sound  speed,  c(r,z)  is  the  acoustic  sound  speed,  and  p(z) 

is  the  density  (which  is  assumed  to  vary  only  with  depth  at  any  given  range).  It  is  within 

c(r,  z)  and  p{z)  that  all  features  of  the  environment  are  represented. 


If  we  now  introduce  the  so-called  PE  field  function  according  to 


'P(r,z)  =  -^j^>(r,z)e 


(2) 


and  assume  that  the  solution  is  dominated  by  the  forward  propagating  field,  it  can  be 
shown  [9]  that  the  2-D  Helmholtz  equation  can  be  reduced  to  a  parabolic  equation  of  the 
form 


dr 


iK(Q-\)'V  , 


(3) 


where  the  square-root  operator,  Q  ,  is  defined  as 


By  defining  the  operators 


Q  = 


2  1  d 


lA 

pdz) 


(4) 


5 


(5) 


s  =  n^  -\ 

^  kl  dz^ 


y. 


1  dp  d 
kip  dz  dz 


the  square -root  operator  may  be  written  as 


Q  =  ^\  +  £  +  p  +  y  . 


(6) 


Solutions  to  the  parabolie  equation  may  then  be  obtained  from  a  one-way 
marehing  algorithm  solution  in  the  general  form 


'F(r-l-Ar,z)  =  exp[/A:oAr(g-l)  'F(r,z)  . 


(V) 


In  practice,  it  is  not  possible  to  directly  apply  the  square-root  operator  in  the  marching 
algorithm  defined  above.  Some  form  of  approximation  is  required. 


Alternatively,  we  could  define  the  field  function  as 

which  then  satisfies  a  PE  of  similar  form 


(8) 


dr 


iK[Q-\)^  . 


In  this  case,  the  square-root  operator  is  defined  by 


(9) 


Q  =  An^  + 


J_iL-  /T 

kl  dz^  ^ 


where 


£  =  n^  -  and  =n^  + 


2kl 


I  gy  3 

p  dz^  2 


^  I  dp^ 
p  dz 


(10) 


(11) 
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The  original  approaches  in  MMPE  were  based  on  this  latter  definition,  which  requires  a 
description  of  piz)  in  terms  of  a  mixing  function  to  smooth  out  the  discontinuity  at  the 
bottom  interface  so  that  the  derivatives  may  be  computed.  The  solution  to  this  PE  is  also 
obtained  from  a  marching  algorithm  of  the  form  given  in  Equation  (7). 

B,  APPROXIMATION  OF  THE  SQUARE-ROOT  OPERATOR 

The  Thomson-Chapman  WAPE  operator  splitting  [8]  produces  an  approximation 
to  the  square-root  operator  in  Equation  (6),  defined  by  [10] 

Ql  =  +  £  +^J^  +  P  +  ^J^  +  /~'^  ■  (12) 

By  using  this  definition,  the  propagator  in  the  marching  algorithm  becomes 

exp [zl^Ar (gj'  - 1)]  =  exp  iS(^^/\  +  £  +  //+7r+7-3) 

«  expj^zA^^/r+y -l)  exp  +  p  -  ij  expj^zA^VTT^-lj  (13) 

=  exp  +  /  -\j  exp  +  -\j  exp[z(5'(n-l)] 

where  S  =  k^/Sr ,  and  the  approximation  in  the  line  is  due  to  the  neglect  of  any  non¬ 
commutation  of  the  operators  s ,  p ,  and  y .  As  noted  by  Yevick  and  Thomson  [10],  the 
order  in  which  the  operations  are  implemented  is  important,  and  the  term  with  y  should 
be  applied  at  the  end  of  the  range  step.  Note  that  the  operator  currently  used  in  MMPE 
corresponds  to 

0i=Vl  +  f  +  p  —  \  ,  (14) 

and  there  is  no  y  term  of  concern. 

The  Thomson-Chapman  approximation  neglects  the  lowest  order  cross  terms 
between  the  operators.  The  lowest  order  correction  is  [10] 


=  0i'-^[^(A+r)+(A+r)^+A7+rA]  • 
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For  completeness,  note  that  the  eorresponding  lowest  order  eorreetion  to  the  existing 
MMPE  operator  would  be 

o 

o 

In  PE  models  based  on  the  SSE  algorithm,  the  operators  eontaining  6  (or  s)  and 
//  are  treated  as  simply  sealar  operators  in  the  physieal  spaee  (z)  and  wavenumber  spaee 

(A:^),  respeetively.  The  veetor  multiplieation  is  then  a  trivial  task,  whieh  leads  to  the 

great  effieieney  of  the  SSE  approaeh.  However,  this  does  require  that  operators 
eontaining  s  and  fJ.  are  separated  in  the  approximation  of  the  square-root  operator. 

If  we  assume  that  the  y  terms  and  eross  terms  ean  be  added  later  as  eorreetions, 
then  an  intermediate  step  solution  ean  be  defined  as 


T'(r  + Ar,z)  =  FFr  '  expi/A  1 — 1--1  i-xFFr|exp[/(5'(n-l)  'P(r,z)|  .(17) 


Note  that  if  we  replaee  n  by  ,  then  this  is  the  eomplete  range  step  utilized  in  the  existing 
MMPE  eode. 


C.  DENSITY  SMOOTHING  APPROACH 

In  order  to  utilize  the  simpler  approaeh  in  MMPE,  based  upon  Equations  (9),  (10), 
(11),  and  (14),  we  must  first  define  the  mixing  funetion  for  the  density  eontrast.  In  this 
ease,  the  density  profile  is  defined  in  terms  of  a  generalized  funetion  representing  the 
diseontinuity  at  the  water-bottom  interfaee  by 

=  (18) 
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is  the  hyperbolic  tangent  mixing  function  approximation  to  the  Heaviside  step  function. 
The  parameter  g  —  z—z^^,  where  z^^  is  the  bottom  depth,  and  L  is  the  mixing  length  of  this 
smoothing  function. 

Now,  since 

Ikl 


1  ay 

p  dz^ 


^  1  dp^ 
p  dz 


,  we  need  to  define 


oz  dz 


4L 


\2Lj 


(20) 


and 


aV  /  ^ 

=  {Pb-P«)- 


{Pb  Pw) 


dz 


dz^ 


41} 


tanh 


v2Ty 


sech^ 


ylLj 


(21) 


Previous  analysis  by  others  [7],  [9],  [11],  [12]  has  shown  that  good  solutions  result  by 

2  ;i 

choosing  a  mixing  length  of  L  =  —  ~  .  Although  this  is  an  additional  free  parameter 

K  ^ 

that  could  be  evaluated  for  accuracy,  we  shall  not  explore  this  further  in  this  thesis  and 
will  keep  this  mixing  length  defined  throughout  for  those  algorithms  that  utilize  the 
density  smoothing  approach. 


D,  YEVICK  AND  THOMSON  HYBRID  IMPLEMENTATION  TECHNIQUE 

The  hybrid  approach  by  Yevick  and  Thomson  [10]  shows  that  the  new  formalism 
defining  QI  requires  evaluation  of  the  additional  operation 


T'(r  +  Ar,z)  =  exp  iS^^jl  +  y  'P(r  +  Ar,z) 


(22) 


Note  that  this  completes  the  range  step  if  we  only  invoke  the  Q'  approximation,  and  we 

may  replace  The  Q2  approximation  requires  an  additional  step  and  is 

discussed  later. 

Yevick  and  Thomson  [10]  proposed  treatment  of  this  operation  through  a  finite 
difference  approach,  whereby  the  Fade  approximation 
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exp 


l  +  -{l  +  iS)/ 


i*\(i-is)r 


(23) 


generates  the  finite  differenee  equations 


1+1(1-, <s)r 


+  Ar,z)  = 


i+i(i+,5)r 


'F(r  +  Ar,z) 


(24) 


dp 

In  this  ease,  direet  applieation  of  y  is  problematie  sinee  y  depends  on  — ,  and 

dz 

we  do  not  want  to  invoke  the  smoothing  funetion  for  p  .  Yeviek  and  Thomson  [10]  show 
that  an  alternative  is  to  evaluate  a  finite  differenee  approximation  to  the  original 
transverse  operator 


/''=4- 

kl  dz 


p  dz 


(25) 


Then  y  ^  p'  -  p  ,  where  p  is  well-defined  as  an  operator  on  'P  . 

Speoifieally,  using  a  3-point  eentered  differenee  approximation  defined  by 


Jo  , 

h 


f 

/.-/( 

V  2  2  y 


+ 


0(1,’-). 


(26) 


then 


/,'T=4A''1«'5' 

kl  dz 


P 


P  )  [k,^z) 


— ('j'.-'j'.)-— ('j'.-'j'-,) 

Px  P  I 


p 


[Ktsz) 


1  . 

1  1 

1 

— 'J^l- 

— + — 

O 

+ 

1 

1 

Px 

'  Px  Px 

P  X 

2 

\  1  ~l) 

2 

(27) 


where  the  indiees  0,  ±  — ,  ±1  apply  to  the  depth  index  in  question  and  neighboring 
depths  above  and  below. 
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Note  that  Yevick  and  Thomson  assumed  a  eontinuous  funetion  for  p  ,  sueh  that 


A  =^(a+a)  ,  P  1  =|(a  +  Ai)  • 

2  ^  2  ^ 

However,  we  shall  treat  the  density  as  defined  by 


P[^)  =  P.^{Pb-  =  \ 


A  ,z<z,; 

|(a  +  a)  .2  =  2,; 


A 


,Z>Zb- 


Then,  sinee 


kl  dz^  kl  dz 


dz 


kg  Az 


az 


2 


az 


/-A 

2 


1 


(kgAz)' 


(kgAz)' 


-['P,-2'Pg+'P_j 


(28) 


(29) 


(30) 


we  obtain 


fV  =  {p'-p)^> 


a  a 

a 

a 

1  1 

a  a 

Po  2 

'I^l- 

Po  1  Po 

-2 

o 

+ 

Po  2 

1 

:  Pi  P  i 

1 

V  2  J 

^2  2 

L  2  J 

•  (31) 


This  agrees  with  Yeviek  and  Thomson  when  their  notation  is  employed,  defined  as 


A  _  A  _  p 


P^ 
P  1 


Po 


1 

2 


(Po+P-l) 


=  p- 


A__i_  _A_ 

Pi  Pi 

2  2 


=  2Po 


-  + 


(Pl+Po)  (Po+P-l) 


=  P++P- 


(32) 


11 


For  our  discontinuous  density,  only  integer  index  values  of  p  have  any  meaning.  So  we 
shall  define  the  notation  /Jq  =  /?  (zg ) ,  Pj  =  p^  =  p{z>  ,  and  p  j  =  p_^  =  p[z  <  z^) . 


Then, 


f¥  = 


(k^Az) 


( 

\ 

(  \ 

~ 

Po  _  ^ 

^1- 

Po  1  Po 

-2 

'i'o  + 

Po  _  Y 

'i'-i 

_vf’+  y 

kP^  P- 

) 

kP-  j 

(33) 


We  ean  then  use  a  tridiagonal  matrix  solution  to  the  expression 


where  A  has  elements  defined  by 


1 


W , 


i+i(i-«)r 


(34) 


and  B  has  elements  defined  by 


i+-(i-,<s)r 


.  Note  that  A  and  B  are  just  identity  matriees  exeept  in  the  vieinity  of 
z  ~  Zi^  ,  where  y  is  non-zero. 


If  we  are  only  interested  in  the  eorreetion  from  the  definition  of  ,  the  solution 

is  eomplete.  However,  Yeviek  and  Thomson  also  go  on  to  deseribe  how  this  ean  be 
further  refined  through  applieation  of  another  finite  differenee  approximation  to  the 

additional  term  in  exp|-^[£(//  +  y)  +  (//  +  y)f +  //y  +  y//]| .  Sinee  p  +  y^p’ , 


we  may  write  this  as 


1 


r  ip  1  l---id{sp'  +  p's  +  py  +  yp) 

expj-— [£•//' +  //'^ +  ///  +  ///]  U — - 

'  1  + —iS[sp' +  p's  +  py  ^-yp') 


(35) 


Note  that  for  y  -  0,  this  eorreetion  simplifies  to 


1 


f  ,1 

expj-— (£•//  + —  1 -  ’ 

^  l  +  —iA(sp  +  ps) 


(36) 


12 


which  can  be  applied  as  a  eorrection  to  the  existing  MMPE  approaeh  based  on  Q2  if  we 
replaee  s^s. 

Yeviek  and  Thomson  drop  the  j^y  and  yju  terms,  arguing  that  these  are  small, 
higher-order  eorrections,  and  instead  utilize  the  eorrection  [10] 


exp 


-y(£’//'  + A) 


l-^i5[sju'  +  ju's) 


(37) 


This  produees  the  additional  eorreetion  for 
'F  (r-l- Ar,z)  =  exp 


Ar,z) 


(38) 


(39) 


where  we  again  invoke  the  matrix  approaeh 

and  elements  of  A  are  defined  by  l+—i5{^£ju'  +  ju'£)  while  the  elements  of  B  are 
defined  by  \——i5{^£ju'  +  ju'£^. 

The  evaluation  of  f/zT  follows  Eq.  (27), 


£ju'^ 


(koAz) 


1 

1  1 

1 

— + — 

'E  +  'E 

^  0  ^  ^  -1 

p\. 

1  Pi  p  1  : 

Pj. 

2 

L  2  2 ) 

2 

(40) 
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However,  for  we  must  consider 


A  A 

"i 

A  5 

■j_r 

kl  dz 

p  dz 

Az  dz 

Po 


(^qAz) 


V  2  2  2  1  ) 


Px  P  1 


(41) 


Po 


(^qAz) 


1  ... 

1  1 

1 

—  ^I'Pi- 

— + — 

^o'Po  + 

Px 

'  Px  Px 

P  X 

2 

\  1  2  J 

2 

Note  that  if  we  are  applying  only  this  correction  to  MMPE,  then  the  correction 
would  result  from  the  application  of  Equation  (36)  to  H'E(r  +  Ar,z)  =  5'E'(r  +  Ar,z) . 
Here, 


(42) 


and 


(43) 


As  before,  we  would  need  to  replace  s^s  if  we  are  utilizing  the  density  smoothing 
approach. 


E,  TAPPERT  SMOOTHING  FUNCTION  APPROACH 


In  1991,  Tappert  [13]  suggested  an  alternative  to  the  effective  index  term  in  the 
smoothing  function  as 


~2  2 
n  =n 


a  d^H{g) 
kl  dz^ 


l  +  (Pw/P(,)2 


(44) 


This  method  can  utilize  the  existing  MMPE  code  by  simply  replacing  the  effective  index 
of  refraction  term  with  Equation  (44).  We  shall  utilize  the  same  definition  for  the  mixing 
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length,  L  =  — ,  although  Tappert  believed  this  should  improve  the  model  as  both  Az  and 

K 

L  get  smaller.  We  shall  not  investigate  the  effeet  of  redueing  the  mixing  length  with  this 
approaeh  further  in  this  thesis. 


15 


THIS  PAGE  INTENTIONALLY  LEET  BLANK 


16 


III.  RESULTS 


A,  MMPE  USER  INTERFACE 

The  MMPE  input  parameters  give  the  user  the  ability  to  manually  define  the  grid 
size.  This  is  accomplished  by  changing  the  “pefiles.inp”  input  file  in  the  local  directory. 
The  depth  grid  size  is  defined  indirectly  by  changing  nz  in  “pefiles.inp”  to  the  desired 
length  of  the  EFT  size  to  be  used  in  the  calculation.  This  value  determines  the  number  of 
depth  samples  at  each  range  step  according  to  Az  =  Iz^^jnz ,  where  z^^^  is  the  maximum 

computational  depth.  The  factor  2  accounts  for  the  image  solution  method  employed  in 
SSF  algorithms.  Similarly,  the  range  step  size  can  be  manually  defined  by  changing  rfact. 
This  defines  the  range  step  by  multiplying  Az  by  rfact.  Therefore,  an  rfact  =  1  makes  a 
uniform  grid  in  which  the  depth  grid  size  and  the  range  step  size  are  equal.  An  rfact  >1 
makes  the  range  step  size  larger  than  the  depth  grid  size  (and  vice  versa).  If  Az  and  rfact 
are  set  to  0.0,  the  model  will  automatically  define  Az  =  /I/IO  and  Ar  =  /I ,  where  X  is  the 
acoustic  wavelength  be  computed. 

B,  NO  DENSITY  CASE 

The  original  MMPE  and  both  alternate  options  were  used  to  model  the 
transmission  loss  in  an  ocean  environment  with  no  density  discontinuity  at  the  bottom 
interface.  In  each  case,  the  MMPE  model  was  allowed  to  automatically  determine  the 
grid  size  based  on  the  default  values  of  Az  and  Ar,  i.e.,  Az  =  /l/10  and  Ar  =  /l.  The 
solution  for  the  alternate  implementation  methods  were  confirmed  to  be  identical  to  the 
solution  obtained  in  the  original  MMPE  model,  since  there  is  no  density  discontinuity 
calculation  in  this  example.  All  solutions  closely  match  the  benchmark  solution,  as 
observed  in  Figure  2. 

Since  both  alternate  implementation  options  exactly  match  the  original  MMPE 
solution  in  this  environment,  only  the  original  MMPE  model  was  used  to  find  the 
optimum  Ar  in  the  case  without  a  density  contrast.  The  Ar  variability  study  was 
conducted  by  varying  rfact  while  holding  nz  constant. 
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No  density  -  accuracy 


30 - 1 - 1 - 1 - 1 - 1 - 1 - r 


100 


110' - ' - ' - ' - ' - ' - ' - ' - ' - ' - 

0  2  4  6  8  10  12  14  16  18  20 

Range  (km) 

Figure  2.  Original  MMPE  vs.  Couple,  no  density  discontinuity. 

a.  Results  for  nz  =  512 

Figures  3-6  show  the  results  for  a  fixed  value  of  nz  =  512,  corresponding  to  a  Az 
of  approximately  5  m,  or  Az«/l/3,  with  rfact  values  of  50,  20,  10,  5,  1,  0.2,  0.1,  and 
0.05.  These  graphs  show  two  notable  features:  that  the  model  converges  to  a  stable 
solution  and  that  the  stable  solution  tends  to  underestimate  the  transmission  loss  of 
the  benchmark  solution  and  lacks  structure.  Focusing  on  the  transmission  loss  from  15- 
20  km,  as  depicted  in  Figure  7,  allows  us  to  better  see  the  model’s  convergence.  For  the 
nz  =  512  case,  the  solution  tends  to  converge  at  rfact  values  smaller  than  1. 

Since  an  rfact  =  5  overestimated  transmission  loss  when  compared  with  the 
benchmark  solution,  and  rfact  =  1  underestimated  the  transmission  loss,  there  must  be  an 
optimum  rfact  between  these  values  which  provides  the  most  accurate  solution  when 
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compared  to  the  benchmark.  The  optimum  value  for  the  nz  =  512  case  is  determined  to  be 
2.5,  as  shown  in  Figure  8.  This  corresponds  to  a  Ar  of  approximately  12  m,  or  Ar  ~  /f . 


MMPE  -  No  density  -  nz  =  512 


Figure  3.  Original  MMPE,  no  density,  nz  =  512,  rfact  =  50,  20. 
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TL  (dB  re  1  m)  for  Pressure  ^rj  TL  (dB  re  1  m)  for  Pressure 


MMPE  -  No  density  -  nz  =  512 


tgure  4. 


Original  MMPE,  no  density,  nz  =  512,  rfact  =  50,  20. 


MMPE  -  No  density  -  nz  =  512 


I  I 


■  Couple 

■  tfact=1 
tfact=.2 


8  10  12  14  16  18  20 
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Figure  5.  Original  MMPE,  no  density,  nz  =  512,  rfact  =  1,  0.2. 
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TL  (dB  re  1m)  for  Pressure 
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MMPE  -  No  density  -  nz  =  512 


Couple 


100 
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0  2  4  6  8  10  12  14  16  18  20 

Range  (km) 

Figure  6.  Original  MMPE,  no  density,  nz  =  512,  rfact  =0.1,  0.05. 


MMPE  -  No  density  -  nz  =  512 


Figure  7.  Solution  eonvergenee  for  nz  =  5 12. 


21 


MMPE  -  No  density  -  nz  =  512 


Figure  8.  Closest  approximation  to  benehmark  solution,  nz  =  5 12. 
b.  Results  for  nz  =  1 024 

The  proeess  was  repeated  as  Azwas  decreased  with  similar  results.  Here,  the 
depth  mesh  is  approximately  10  m,  or  Az  ~  /1/6 .  The  solution  again  converges  at  a  value 
that  slightly  underestimates  transmission  loss  as  observed  in  Figures  9-13.  In  this  case, 
the  optimum  rfact  is  determined  to  be  5  as  shown  in  Figure  14.  This  again  corresponds  to 
a  Ar  of  approximately  12  m,  or  Ar  ~  /I . 
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MMPE-Nodensity-nz-  1024 


Figure  9.  Original  MMPE,  no  density,  nz  =  1024,  rfact  =  50,  20. 


MMPE-No(jen»i(y-nz-  1024 


10  12  14  18  18  20 


R«nge(km) 


Figure  10.  Original  MMPE,  no  density,  nz  =  1024,  rfact  =  10,  5. 
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MMPE-Nodenjity-nz-  1024 


Figure  11.  Original  MMPE,  no  density,  nz  =  1024,  rfact  =  1,  0.2. 


MMPE-No(len$ily-n2»  1024 


Figure  12.  Original  MMPE,  no  density,  nz  =  1024,  rfact  =  0.1,  0.05. 
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MMPE  -  No  density  -  nz  =  1024 


Figure  13.  Solution  convergence  for  nz  =  1024. 


MMPE  -  No  density  -  nz  =  1024 


Figure  14.  Closest  approximation  to  benchmark  solution,  nz  =  1024. 
c.  Results  for  nz>1024 

The  model  was  also  evaluated  at  nz  values  of  2048,  4096,  8192,  and  16,384.  The 
results  for  smaller  Azwere  similar  to  those  obtained  for  nz  =  512  and  1024.  In  all  cases, 
the  solution  converges  at  a  value  slightly  lower  than  the  benchmark  solution  and  the 
optimum  solution  occurs  when  Aris  approximately  12  m,  or  Ar~/l.  See  the  Appendix 
for  full  results. 
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c. 


MODEL  RESULTS  IN  PRESENCE  OF  DENSITY  DISCONTINUITY 


All  three  implementation  options  for  the  MMPE  model  were  employed  in  an 
environment  with  a  density  contrast  at  the  bottom  interface.  The  range  step  size  was  held 
constant  at  12  m  while  the  depth  grid  size  was  reduced  for  each  implementation  method. 
Range  step  correction  factors  slightly  higher  and  lower  were  also  tested  against  a  varying 
depth  grid  size. 

1,  Original  MMPE 

With  a  fixed  Ar=12  m,  there  are  negligible  differences  in  the  transmission  loss 
model  as  observed  in  Figure  15.  In  this  case,  an  nz  value  of  512  is  not  only  the  most 
efficient  calculation,  but  also  provides  a  transmission  loss  solution  that  is  not 
significantly  different  than  those  obtained  from  a  finer  depth  grid.  As  noted,  for 
nz  =  512 ,  the  depth  mesh  is  approximately  Az  ~  5  m  ~  /1/3 . 


Original  MMPE  -  with  density  -  deltar  =  12  m 


Figure  15.  Original  MMPF,  Ar=  12  m,  various  nz. 


2,  MMPE  with  Tappert  Smoothing  Function 

Figure  16  shows  similar  results  when  the  Tappert  smoothing  function  is 
implemented  into  the  MMPF  model.  When  compared  to  the  original  MMPF,  the 
differences  in  the  model  are  even  less  significant  as  Az  is  decreased. 
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MMPE  with  correction  -  with  density  -  deltar  =  1 2  m 


Figure  16.  MMPE  with  Tappert  smoothing  function,  Ar=  12  m,  various  nz. 

3,  Yevick  and  Thomson  Hybrid  Implementation 

Yevick  and  Thomson’s  SSF/FD  hybrid  implementation  provides  for  some 
interesting  results.  At  the  largest  Az  (nz  =  512),  the  model  is  a  poor  comparison  to  the 
benchmark  solution  at  range.  The  solution  underestimates  the  benchmark  solution  at  all 
values.  Furthermore  the  shape  of  the  graph  only  vaguely  represents  the  benchmark 
solution.  Decreasing  the  depth  grid  size  by  a  factor  of  2  results  in  an  even  worse  solution 
as  seen  in  Figure  17.  Flowever,  decreasing  Az  beyond  this  point  improved  the  solution 
dramatically.  At  nz  =  2048,  the  solution  begins  approaching  the  benchmark  solution,  both 
in  overall  shape  and  magnitude  of  transmission  loss.  The  hybrid  model  most  closely 
approximates  the  benchmark  solution  at  nz  =  8192,  corresponding  to  Az  «  0.3  m  -  2/50  . 
Decreasing  Az  further  results  in  smaller  changes  to  the  solution  as  observed  in  Figure  18. 
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SSF/FD  hybrid  -  with  density  -  deltar  =  1 2  m 


Figure  17.  Flybrid  implementation,  Ar=  12  m,  various  nz  (coarse  Az). 


SSF/FD  hybrid  -  with  density  -  deltar  =  1 2  m 


Figure  18.  Hybrid  implementation,  Ar=  12  m,  various  nz  (fine  Az). 

4,  Overall  Comparison 

Solutions  for  the  original  MMPE,  MMPE  with  Tappert  smoothing  function,  and 
Yevick  and  Thomson’s  hybrid  implementation  are  compared  in  Eigure  19.  An  nz  value  of 
512(Az~/l/3)is  used  for  both  the  original  version  of  the  MMPE  model  and  the  MMPE 
model  with  the  Tappert  alternative  method,  since  that  value  was  efficient  and  smaller  Az 
values  did  not  significantly  affect  the  solution  accuracy.  An  nz  value  of  8 1 92  ( Az  ~  /1/50  ) 
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is  used  for  the  hybrid  implementation  option  sinee  that  value  provided  for  the  most 
accurate  solution  for  the  hybrid  implementation  technique.  The  figure  shows  there  is 
marginal  difference  between  the  original  MMPE  and  MMPE  with  the  Tappert  alternative 
function.  The  hybrid  implementation  option  provides  for  a  slightly  closer  approximation 
to  the  benchmark  solution  in  amplitude,  however  is  a  less  efficient  solution.  In  all  cases, 
though,  the  phase  error  exists  at  this  long  range.  Thus,  it  does  not  appear  that  this  phase 
error  is  caused  by  inaccuracies  in  the  computational  treatment  of  the  smoothing  function, 
or  the  first  order  correction  in  the  hybrid  method  defined  in  Equation  (22). 


Implementation  Comparison  -  with  density  -  deltar  =  12  m 


Eigure  19.  Implementation  method  comparison. 


5,  Effect  of  Reference  Sound  Speed 

To  this  point,  the  reference  sound  speed  for  all  transmission  loss  solutions  was 
held  constant  at  1500  m/s.  However,  varying  the  reference  sound  speed  has  a  significant 
effect  on  the  phase  of  the  solution.  In  general,  the  magnitudes  of  the  solutions  are 
unchanged,  but  the  phase  error  present  can  be  minimized.  Eigures  20-23  show  the  results 
when  the  reference  sound  speed  is  varied.  The  previously  determined  optimum  depth 
mesh  and  range  step  is  used  for  each  implementation  option.  The  best  solution  appears  to 
occur  at  a  reference  sound  speed  of  1480  m/s.  Eigure  24  shows  the  results  of  all  three 
implementation  options  at  the  optimum  reference  sound  speed. 
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TL  (dB  re  1m)  for  Pressure  TL  (dB  re  1m)  for  Pressure 


Original  MM  PE  -  with  density  -  deltar  =  1 2  m  -  various  c 


Range  (km) 


Figure  20.  Original  MMPE,  Ar  =  12  m,  nz  =  5 12,  various  Cq. 


MM  PE  with  correction  -  with  density  -  deltar  =  1 2  m  -  various  c 


Figure  21.  MMPE  with  Tapper!  smoothing,  Ar  =  12  m,  nz  =  512,  various  Cq. 
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SSF/FD  hybrid  -  with  density  -  deltar  =  1 2  m  -  various  c 


Range  (km) 


Figure  22.  Hybrid  implementation,  Ar  =  12  m,  nz  =  8192,  various  Cq. 


Implementation  comparison  -  with  density  -  deltar  =  1 2  m  -  c0=1480 


Figure  23.  Implementation  comparison,  Ar 


12  m,  15-20  km. 
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Implementation  comparison  -  with  density-  deltar  =  12  m  -  c0=1480 
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Range  (km) 


Figure  24.  Implementation  comparison,  Ar  =  12  m,  0-20  km. 


The  model’s  inaccuracies  at  range  appear  primarily  due  to  the  phase  errors. 
Alternate  implementation  options  did  not  significantly  correct  the  phase  error.  Proper 
selection  of  reference  sound  speed  seems  to  correct  this  error.  This  was  previously  noted 
in  a  deep  water  benchmark  problem  defined  in  the  1993  PE  Workshop  [12],  and 
subsequent  work  by  Tappert  [13]  using  his  “co-insensitive”  version  seemed  to  correct  the 
phase  errors.  Yevick  and  Thomson  [10]  were  also  able  to  correct  the  phase  error  found  in 
the  deep  water  PE  Workshop  problem  by  implementing  their  higher-order  correction 
term,  defined  in  Equation  (38),  in  their  hybrid  method.  It  is  quite  possible  that 
implementation  of  that  higher-order  correction  term  would  also  address  the  phase  error 
found  in  the  shallow  water  problem  defined  in  this  thesis.  Although  not  as  efficient  as  a 
basic  SSF  method,  the  proposed  hybrid  method  with  higher-order  correction  might  be 
expected  to  improve  the  long-range  phase  accuracy  of  all  solutions. 
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IV.  SUMMARY 


The  parabolic  equation  method  for  transmission  loss  modeling  is  an  important 
topic  in  the  acoustics  scientific  community.  The  SSF  technique  has  been  used  for  years 
due  to  its  computational  efficiency.  Presently,  SSF  techniques  provide  accurate  solutions 
in  deep  ocean  environments  where  bottom  interactions  are  negligible.  However,  as  the 
focus  shifts  to  shallow  water  acoustic  propagation,  it  becomes  necessary  to  find  a  way  to 
use  the  SSF  algorithm  to  model  transmission  loss.  This  paper  compares  two  different 
implementation  techniques  when  compared  to  a  standard  SSF  model. 

The  basic  SSF  algorithm  employing  the  Thomson-Chapman  WAPE 
approximation  employed  in  the  MMPE  model  was  found  to  be  very  accurate  when  there 
is  no  density  discontinuity  present.  As  the  depth  mesh  grid  is  reduced  in  size,  the  MMPE 
model  converges  to  a  stable  solution.  When  compared  to  the  benchmark  solution,  the 
convergence  in  Ar  slightly  underestimates  transmission  loss  at  all  ranges.  This  is 
presumably  due  to  numerical  noise  when  the  range  step  becomes  too  small.  By 
systematically  reducing  Ar  it  is  shown  that  an  optimal  comparison  exists  at  a  Ar  ~  1  . 
This  optimal  value  for  Ar  remains  constant  even  as  Az  is  decreased  and  a  density 
discontinuity  is  introduced. 

In  the  presence  of  a  density  discontinuity,  both  the  original  MMPE,  and  MMPE 
with  the  alternate  Tappert  approach  to  treating  the  effective  index  of  refraction  function, 
provide  a  fair  model  of  the  magnitude  of  transmission  loss  at  ranges  greater  than  15  km. 
The  largest  tested  Az  of  approximately  /1/4  provided  acceptable  results.  Decreasing  Az 
further  does  not  significantly  change  the  solutions’  accuracy.  Since  a  smaller  Az  results 
in  increased  computational  time,  it  is  recommended  that  a  Az  of  approximately  /1/4  be 
used  for  both  options. 

The  implementation  of  Yevick  and  Thomson’s  hybrid  technique  was  able  to 
provide  a  slightly  more  accurate  comparison  to  the  benchmark  solution  amplitude. 
However,  it  is  necessary  to  make  Az  very  small  in  order  to  obtain  such  comparable 
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results.  This  comes  at  the  cost  of  increased  computational  complexity,  and  should  be 
balanced  against  the  need  to  obtain  the  most  accurate  result. 

All  models  were  found  to  suffer  from  a  phase  error,  which  accumulates  in  range 
when  a  reference  sound  speed  of  1500  m/s  is  used.  This  phase  error  can  be  minimized  by 
lowering  co  to  1480  m/s.  This  is  reminiscent  of  previous  phase  errors  found  with  the 
Thomson-Chapman  WAPE  approximation.  It  is  quite  possible  that  these  phase  errors 
may  be  significantly  reduced  by  employing  the  higher-order  correction  of  Yevick  and 
Thomson,  or  the  Co-insensitive  version  of  the  SSF  algorithm  defined  by  Tappert.  Future 
work  could  investigate  the  relative  accuracy  of  each  of  these  methods. 
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APPENDIX.  MMPE  CONVERGENCE  GRAPHS 


MMPE  -No  Oenjity-nz  •  2048 
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Figure  25.  Original  MMPE,  no  density,  nz  =  2048,  rfact  =  50,  20. 
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Figure  26.  Original  MMPE,  no  density,  nz  =  2048,  rfact  =  10,  5. 
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MMPE  -No  denoily-nz  -  2048 
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Figure  27.  Original  MMPE,  no  density,  nz  =  2048,  rfact  =  1,  0.2. 


MMPE -No  density-nz  ■  2048 


Figure  28.  Original  MMPE,  no  density,  nz  =  2048,  rfact  =  0.1,  0.05. 
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TL  (dB  re  1m)  for  Pressure  TL  (dB  re  1m)  for  Pressure 


MMPE  -  No  density  -  nz  =  2048 


Range  (km) 


Figure  29.  Solution  convergence  for  nz  =  2048. 


MMPE  -  No  density  -  nz  =  2048 


Figure  30.  Closest  approximation  to  benchmark  solution,  nz  =  2048. 
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MMPE  -  No  density  -  nz  =  4096 
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Figure  3 1 .  Original  MMPE,  no  density,  nz  =  4096,  rfact  =  100,  50. 
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Figure  32.  Original  MMPE,  no  density,  nz  =  4096,  rfact  =  20,  10. 
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TL  (dB  re  1m)  for  Pressure 
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Figure  33.  Original  MMPE,  no  density,  nz  =  4096,  rfact  =  5,  1. 


MMPE  -  No  density  -  nz  =  4096 


Figure  34.  Solution  convergence  for  nz  =  4096. 
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TL  (dB  re  1  m)  for  Pressure 


MMPE  -  No  density  -  nz  =  4096 


Figure  35.  Closest  approximation  to  benchmark  solution,  nz  =  4096. 
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Figure  36.  Original  MMPE,  no  density,  nz  =  8192,  rfact  =  200,  100. 
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Figure  37.  Original  MMPE,  no  density,  nz  =  8192,  rfact  =  50,  20. 
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Figure  38.  Original  MMPE,  no  density,  nz  =  8192,  rfact  =  10,  5. 
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TL  (dB  re  1m)  for  Pressure  TL  (dB  re  1m)  for  Pressure 
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Figure  39.  Solution  convergence  for  nz  =  8192. 


MMPE  -  No  density  -  nz  =  8192 
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Figure  40.  Closest  approximation  to  benchmark  solution,  nz  =  8192. 
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Figure  41.  Original  MMPE,  no  density,  nz  =  16,384,  rfact  =  800,  400. 
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Figure  42.  Original  MMPE,  no  density,  nz  =  16,384,  rfact  =  200,  100. 
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Figure  43.  Original  MMPE,  no  density,  nz  =  16,384,  rfact  =  50,  20. 
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Figure  44.  Solution  convergence  for  nz  =  16,384. 
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TL  (dB  re  1  m)  for  Pressure 


MMPE  -  No  density  -  nz  =  1 6384 


Figure  45.  Closest  approximation  to  benchmark  solution,  nz  =  16,384. 
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